/*
 * DoubleSparseMatrix.cs
 * 
 * Copyright (c) 2003-2004, dnAnalytics Project. All rights reserved.
*/

using System;
using System.Collections.Generic;
using System.Diagnostics.CodeAnalysis;
using dnAnalytics.Exceptions;
using dnAnalytics.Resources;

namespace dnAnalytics.LinearAlgebra
{
    /// <summary>
    /// Defines a matrix class that only stores non-zero values.
    /// </summary>
    /// <remarks>
    /// In the current implementation the <b>SparseMatrix</b> uses a
    /// compressed row storage system.
    /// </remarks>
    internal sealed class SparseMatrix : Matrix
    {
        // UPGRADE NOTE: Could change the current storage scheme from
        // Compressed Row Storage (currently implemented) to a special
        // row storage scheme which uses individual arrays for each row
        // There would be one object array for the row 'list'
        // Each item in the row array would contain either null (no values in that row)
        // or an array of structs (the struct contains the value and the 
        // column index) which acts as a column array. Only non zero values are stored
        // This way rowwise access may be a bit quicker and as the columns 
        // will in general be smaller than a full column they might also
        // fit in the CPU cache better.

        // Can also use sorted lists for the sub columns.

        /// <summary>
        /// The array containing the actual non-zero values. Stored in a compressed
        /// row fashion.
        /// </summary>
        internal double[] mValues = new double[0];

        /// <summary>
        /// An array containing the column indices of the non-zero values.
        /// </summary>
        internal int[] mColumnIndices = new int[0];

        /// <summary>
        /// The array containing the row indices of the existing rows. 
        /// </summary>
        internal int[] mRowIndices = new int[0];

        /// <summary>
        /// The total number of non-zero values in the matrix.
        /// </summary>
        /// <remarks>
        /// This variable keeps track of the number of non-zero values that are stored
        /// in the <b>mValues</b> and <b>mColumnIndices</b> array's. Because of this
        /// <b>mValueCount</b> also indicates the 'length' of the array's. Note that the
        /// actual arrays may be longer (but never shorter) than indicated by
        /// <b>mValueCount</b>. If the actual arrays are longer than indicated reading
        /// beyond the number of values (i.e. <b>mValueCount</b>) will produce garbage
        /// data.
        /// </remarks>
        private int mValueCount;

        /// <summary>
        /// Creates a new instance of the <see cref="SparseMatrix"/> class.
        /// </summary>
        /// <param name="dimension">The number of rows and columns in the matrix.</param>
        public SparseMatrix(int dimension)
            : this(dimension, dimension)
        {
            // Redirects to SparseMatrix(int rows, int columns)
        }

        /// <summary>
        /// Creates a new instance of the <see cref="SparseMatrix"/> class.
        /// </summary>
        /// <param name="rows">The number of rows in the matrix.</param>
        /// <param name="columns">The number of columns in the matrix.</param>
        public SparseMatrix(int rows, int columns)
            : base(rows, columns)
        {
            mRowIndices = new int[rows];
        }

        /// <summary>
        /// Creates a new instance of the <see cref="SparseMatrix"/> class.
        /// </summary>
        /// <remarks>
        /// When creating a new <see cref="SparseMatrix"/> based on another matrix
        /// one should consider carefully how many non-zero values there are in the
        /// <paramref name="source"/> matrix. The useage of a sparse matrix only makes
        /// sense if the matrix is actually 'sparse', i.e. contains relatively few 
        /// non-zero values.
        /// </remarks>
        /// <param name="source">
        /// The <see cref="Matrix"/> from which the data will
        /// be copied to the new instance.
        /// </param>
        public SparseMatrix(Matrix source)
            : base(source.Rows, source.Columns)
        {
            mRowIndices = new int[Rows];

            for (int i = 0; i < source.Rows; i++)
            {
                for (int j = 0; j < source.Columns; j++)
                {
                    InsertValue(i, j, source[i, j]);
                }
            }
        }

        /// <summary>
        /// Creates a new instance of the <see cref="SparseMatrix"/> class.
        /// </summary>
        /// <param name="source">
        /// The <see cref="SparseMatrix"/> from which the data will
        /// be copied to the new instance.
        /// </param>
        public SparseMatrix(SparseMatrix source)
            : base(source.Rows, source.Columns)
        {
            mRowIndices = new int[source.mRowIndices.Length];

            // Copy the data
            source.CopyToSparseMatrix(this);
        }

        /// <summary>
        /// Creates a new instance of the <see cref="SparseMatrix"/> class.
        /// </summary>
        /// <remarks>
        /// When creating a new <see cref="SparseMatrix"/> based on an array of doubles
        /// one should consider carefully how many non-zero values there are in the
        /// <paramref name="source"/>. The useage of a sparse matrix only makes
        /// sense if the matrix is actually 'sparse', i.e. contains relatively few 
        /// non-zero values.
        /// </remarks>
        /// <param name="source">The array of doubles from which the data will be
        /// copied to the new instance.</param>
        public SparseMatrix(double[,] source)
            : base(source.GetLength(0), source.GetLength(1))
        {
            mRowIndices = new int[Rows];

            for (int i = 0; i < source.GetLength(0); i++)
            {
                for (int j = 0; j < source.GetLength(1); j++)
                {
                    InsertValue(i, j, source[i, j]);
                }
            }
        }

        // UPGRADE NOTE: Could potentially store last accessed item (the index numbers
        // because most of the time matrices are accessed sequentially which means
        // that we can figure out what item to pick without having to actually
        // search for it.

        /// <summary>
        /// Finds either the item indicated by the given indices (returns <b>true</b>) or the 
        /// insert position (returns <b>false</b>) in the values array. 
        /// </summary>
        /// <param name="rowIndex"></param>
        /// <param name="columnIndex"></param>
        /// <param name="itemIndex"></param>
        /// <returns></returns>
        private bool FindItem(int rowIndex, int columnIndex, out int itemIndex)
        {
            if ((rowIndex < 0) || (rowIndex >= Rows) || (columnIndex < 0)
                || (columnIndex >= Columns))
            {
                itemIndex = -1;
                return false;
            }

            // First find the start and end indices for the row
            int startIndex;
            int endIndex;
            FindRowStartAndEnd(rowIndex, out startIndex, out endIndex);

            // Check if startIndex and endIndex are equal. If so then
            // the item does not exist
            if (startIndex == endIndex)
            {
                itemIndex = endIndex;
                return false;
            }

            // Now subtract 1 from the endIndex because the current endIndex
            // indicates were the next row starts. 
            endIndex -= 1;

            // If the startIndex is only one bigger than the endIndex
            // do a seperate check because the binary search loop will not 
            // be entered.
            // Check if one of the boundaries is the correct index
            if (startIndex >= endIndex - 1)
            {
                return FindItemForSizeTwoArray(startIndex, endIndex, columnIndex, out itemIndex);
            }

            itemIndex = (startIndex + endIndex) / 2;
            // Start the binary search
            while (startIndex < endIndex - 1)
            {
                if (mColumnIndices[itemIndex] == columnIndex)
                {
                    // Found the value we are looking for so we can exit
                    // here
                    return true;
                }

                // We're still here so we didn't find it. Now see in
                // what direction to search
                if (mColumnIndices[itemIndex] > columnIndex)
                {
                    // value is to big that means that we need to move 
                    // back in the array
                    endIndex = itemIndex;
                }
                else
                {
                    // So the value must be to small
                    // Move forward in the array
                    startIndex = itemIndex;
                }
                itemIndex = (startIndex + endIndex) / 2;

                // Do a quick check to see if the end of the search was
                // reached. This occurs when itemIndex is equal to one
                // of the two boundaries. In this case this should always
                // be the lower end because integer division trunkates instead
                // of rounding
                if (startIndex >= endIndex - 1)
                {
                    return FindItemForSizeTwoArray(startIndex, endIndex, columnIndex, out itemIndex);
                }
            }
            return false;
        }

        /// <summary>
        /// Locates the start and end of the row with the specified index. 
        /// </summary>
        /// <param name="rowIndex">The index of the row that must be located.</param>
        /// <param name="startIndex">The index of the first element in the <see cref="mValues"/> array that belongs
        /// to the given row.</param>
        /// <param name="endIndex">The index of the first element in the <see cref="mValues"/> array that belongs
        /// to the given row + 1.</param>
        private void FindRowStartAndEnd(int rowIndex, out int startIndex, out int endIndex)
        {
            // First find the start and end indices for the row
            startIndex = mRowIndices[rowIndex];
            if (rowIndex < mRowIndices.Length - 1)
            {
                endIndex = mRowIndices[rowIndex + 1];
            }
            else
            {
                // We are looking for a value in the last row so
                // the end of the search is reached as we get to the
                // last valid item in the column array
                endIndex = mValueCount;
            }
        }

        private bool FindItemForSizeTwoArray(int startIndex, int endIndex, int columnIndex, out int itemIndex)
        {
            if ((startIndex > -1) && (startIndex < mColumnIndices.Length)
                &&
                (mColumnIndices[startIndex] == columnIndex))
            {
                itemIndex = startIndex;
                return true;
            }

            if ((endIndex > -1) && (endIndex < mColumnIndices.Length)
                &&
                (mColumnIndices[endIndex] == columnIndex))
            {
                itemIndex = endIndex;
                return true;
            }

            // neither one of the boundaries is actually correct
            // Check where eventual new values should be placed
            if (mColumnIndices[startIndex] > columnIndex)
            {
                itemIndex = startIndex;
            }
            else
            {
                if (mColumnIndices[endIndex] > columnIndex)
                {
                    itemIndex = endIndex;
                }
                else
                {
                    itemIndex = endIndex + 1;
                }
            }
            return false;
        }

        private void InsertValue(int rowIndex, int columnIndex, double value)
        {
            // Check if the value is zero. If it is then do nothing 
            if (System.Math.Abs(value)
                > Double.Epsilon)
            {
                // We assume that the item is not present in the array
                // First find the spot where the value must be inserted
                int index;
                FindItem(rowIndex, columnIndex, out index);

                // Check if the storage needs to be increased
                if (mValueCount == mValues.Length)
                {
                    // Value array is completely full so we increase the size
                    // This also grows the columIndex array
                    GrowStorage();
                }
                // Move all values (with an position larger than index) in the value array 
                // to the next position
                // move all values (with an position larger than index) in the columIndices 
                // array to the next position
                for (int i = mValueCount - 1; i > index - 1; i--)
                {
                    mValues[i + 1] = mValues[i];
                    mColumnIndices[i + 1] = mColumnIndices[i];
                }

                // Add the value and the column index
                mValues[index] = value;
                mColumnIndices[index] = columnIndex;

                // increase the number of non-zero numbers by one
                mValueCount += 1;

                // add 1 to all the row indices for rows bigger than rowIndex
                // so that they point to the correct part of the value array 
                // again.
                for (int i = rowIndex + 1; i < mRowIndices.Length; i++)
                {
                    mRowIndices[i] += 1;
                }
            }
        }

        private void DeleteItem(int rowIndex, int columnIndex)
        {
            // We have an item that needs to be removed from the
            // array. First check if it is actually stored in the array
            int index;
            if (FindItem(rowIndex, columnIndex, out index))
            {
                // Move all values (with an position larger than index) in the value array 
                // to the previous position
                // move all values (with an position larger than index) in the columIndices 
                // array to the previous position
                for (int i = index + 1; i < mValueCount; i++)
                {
                    mValues[i - 1] = mValues[i];
                    mColumnIndices[i - 1] = mColumnIndices[i];
                }

                for (int i = rowIndex + 1; i < mRowIndices.Length; i++)
                {
                    mRowIndices[i] -= 1;
                }
                mValueCount -= 1;
            }
        }

        private void GrowStorage()
        {
            // Grow the valueArray by X%
            // Grow the columnIndex array by X%
            // Both arrays should always be equal in length!!!!!
            int size = mValues.Length + GrowthSize();
            double[] newValues = new double[size];
            int[] newIndices = new int[size];
            for (int i = 0; i < mValues.Length; i++)
            {
                newValues[i] = mValues[i];
                newIndices[i] = mColumnIndices[i];
            }
            mValues = newValues;
            mColumnIndices = newIndices;
        }

        private int GrowthSize()
        {
            int delta;
            if (mValues.Length > 1024)
            {
                delta = mValues.Length / 4;
            }
            else
            {
                if (mValues.Length > 256)
                {
                    delta = 512;
                }
                else
                {
                    if (mValues.Length > 64)
                    {
                        delta = 128;
                    }
                    else
                    {
                        delta = 32;
                    }
                }
            }
            return delta;
        }

        protected override Matrix CreateMatrix(int rows, int columns)
        {
            return MatrixBuilder.CreateMatrix(rows, columns, MatrixType.Sparse);
        }

        /// <summary>
        /// Retrieves the value at the given indices. Note: This method is not ranged checked. If the
        /// indices are outside the dimensions of the matrix, it may return an incorrect value.
        /// </summary>
        /// <param name="row">The row index.</param>
        /// <param name="column">The column index.</param>
        /// <returns>The value at the given indices.</returns>
        protected internal override double ValueAt(int row, int column)
        {
            int index;
            if (FindItem(row, column, out index))
            {
                return mValues[index];
            }
            return 0.0;
        }

        /// <summary>
        /// Set the value at the given indices to the given value. Note: This method is not ranged checked. If the
        /// indices are outside the dimensions of the matrix, it may set an incorrect element or throw
        /// an <see cref="IndexOutOfRangeException"/>.
        /// </summary>
        /// <param name="row">The row index.</param>
        /// <param name="column">The column index.</param>
        /// <param name="value">The value to set.</param>
        /// <exception cref="IndexOutOfRangeException">Depending on the implementation, an <see cref="IndexOutOfRangeException"/>
        /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
        protected internal override void ValueAt(int row, int column, double value)
        {
            int index;
            if (FindItem(row, column, out index))
            {
                if (System.Math.Abs(value)
                    < Double.Epsilon)
                {
                    DeleteItem(row, column);
                }
                else
                {
                    mValues[index] = value;
                }
            }
            else
            {
                InsertValue(row, column, value);
            }
        }

        /// <summary>
        /// Returns the number of non zero elements in the matrix.
        /// </summary>
        /// <value>The number of non zero elements.</value>
        public int NonZeros
        {
            get { return mValueCount; }
        }

        /// <summary>Calculates the Frobenius norm of this matrix.</summary>
        /// <returns>The Frobenius norm of this matrix.</returns>
        public override double FrobeniusNorm()
        {
            // norm = sqrt(the sum of the squares of all items)
            double norm = 0;
            for (int i = 0; i < mValueCount; i++)
            {
                norm += mValues[i] * mValues[i];
            }

            return System.Math.Sqrt(norm);
        }

        /// <summary>Calculates the infinity norm of this matrix.</summary>
        /// <returns>The infinity norm of this matrix.</returns>      
        public override double InfinityNorm()
        {
            // The infinity norm is the largest row sum of the matrix
            // This one is simple because all is sorted by rows
            double norm = double.NegativeInfinity;
            for (int i = 0; i < Rows; i++)
            {
                int startIndex;
                int endIndex;
                FindRowStartAndEnd(i, out startIndex, out endIndex);

                double rowSum = 0;
                for (int j = startIndex; j < endIndex; j++)
                {
                    rowSum += System.Math.Abs(mValues[j]);
                }
                norm = System.Math.Max(norm, rowSum);
            }
            return norm;
        }

        /// <summary>
        /// Sets all values to zero.
        /// </summary>
        public override void Clear()
        {
            mValueCount = 0;
            Array.Clear(mRowIndices, 0, mRowIndices.Length);
        }

        /// <summary>
        /// Copies a column into a vector.
        /// </summary>
        /// <param name="index">The column to copy.</param>
        /// <returns>A vector containing the copied elements.</returns>
        public override Vector Column(int index)
        {
            Vector result = VectorBuilder.CreateVector(Rows, VectorType.Sparse);
            Column(index, result);
            return result;
        }

        /// <summary>
        /// Copies a column into to the given vector.
        /// </summary>
        /// <param name="index">The column to copy.</param>
        /// <param name="result">The vector to copy the column into.</param>
        /// <exception cref="ArgumentNullException">If the result matrix is null.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="index"/> is negative,
        /// or greater than or equal to the number of columns.</exception>        
        /// <exception cref="NotConformableException">If the result vector is of the incorrect size.</exception>
        public override void Column(int index, Vector result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (index >= Columns || index < 0)
            {
                throw new ArgumentOutOfRangeException("index", Strings.OutOfRangeException);
            }

            if (result.Count != Rows)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            
            // Clear the vector
            result.Clear();

            for (int i = 0; i < Rows; i++)
            {
                int itemIndex;
                if (!FindItem(i, index, out itemIndex))
                {
                    continue;
                }

                result[i] = mValues[itemIndex];
            }
        }

        /// <summary>
        /// Copies the given <see cref="Range"/> of a column into a vector.
        /// </summary>
        /// <param name="index">The column to copy.</param>
        /// <param name="range">The <see cref="Range"/> of elements to copy.</param>
        /// <returns>A vector containing the copied elements.</returns>
        public override Vector Column(int index, Range range)
        {
            // the actual size of the vector is the number of elements in the 
            // range. This should be the difference between start and end + 1  (<-- off by one thing)
            // and then the number of times the stride fits in that difference.
            int size = range.End - range.Start;
            size = size / range.Stride + 1;

            Vector result = VectorBuilder.CreateVector(size, VectorType.Sparse);
            Column(index, range, result);
            return result;
        }

        /// <summary>
        /// Copies the given <see cref="Range"/> of a column into the given vector.
        /// </summary>
        /// <param name="index">The column to copy.</param>
        /// <param name="range">The <see cref="Range"/> of elements to copy from the requested column.</param>
        /// <param name="result">The <see cref="Vector"/> to copy the column into.</param>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="range"/> specifies a range
        /// outside the column's range. </exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <b>null</b>.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="index"/> is negative,
        /// or greater than or equal to the number of columns.</exception>        
        /// <exception cref="NotConformableException">If <b>result.Count != (range.End - range.Start) / range.Stride + 1</b>.</exception>
        public override void Column(int index, Range range, Vector result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }
            if (index >= Columns || index < 0)
            {
                throw new ArgumentOutOfRangeException("index", Strings.OutOfRangeException);
            }
            if (range.Start >= Rows || range.End >= Rows)
            {
                throw new ArgumentOutOfRangeException("range", Strings.OutOfRangeException);
            }
            
            // Clear the result vector
            result.Clear();

            for (int i = range.Start, k = 0; i <= range.End; i += range.Stride, k++)
            {
                int itemIndex;
                if (!FindItem(i, index, out itemIndex))
                {
                    continue;
                }

                result[k] = mValues[itemIndex];
            }
        }

        /// <summary>
        /// Copies a row into a vector.
        /// </summary>
        /// <param name="index">The row to copy.</param>
        /// <returns>A vector containing the copied elements.</returns>
        public override Vector Row(int index)
        {
            Vector result = VectorBuilder.CreateVector(Columns, VectorType.Sparse);
            Row(index, result);
            return result;
        }

        /// <summary>
        /// Copies a row into to the given vector.
        /// </summary>
        /// <param name="index">The row to copy.</param>
        /// <param name="result">The vector to copy the row into.</param>
        /// <exception cref="ArgumentNullException">If the result matrix is null.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="index"/> is negative,
        /// or greater than or equal to the number of rows.</exception>        
        /// <exception cref="NotConformableException">If the result vector is of the incorrect size.</exception>
        public override void Row(int index, Vector result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (index >= Rows
                || index < 0)
            {
                throw new ArgumentOutOfRangeException("index", Strings.OutOfRangeException);
            }

            if (result.Count != Columns)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            
            // Clear the result vector
            result.Clear();

            // The data is stored in compressed row storage so all we have to do is find the
            // begin and end of the row and we can quickly iterate over that.
            int begin;
            int end;
            FindRowStartAndEnd(index, out begin, out end);

            for (int i = begin; i < end; i++)
            {
                int insertPosition = mColumnIndices[i];
                result[insertPosition] = mValues[i];
            }
        }

        /// <summary>
        /// Copies the given <see cref="Range"/> of a row into a vector.
        /// </summary>
        /// <param name="index">the column to copy.</param>
        /// <param name="range">the <see cref="Range"/> of elements to copy.</param>
        /// <returns>a vector containing the copied elements.</returns>
        public override Vector Row(int index, Range range)
        {
            // the actual size of the vector is the number of elements in the 
            // range. This should be the difference between start and end + 1 (<-- off by one thing)
            // and then the number of times the stride fits in that difference.
            int size = range.End - range.Start;
            size = size / range.Stride + 1;

            Vector result = VectorBuilder.CreateVector(size, VectorType.Sparse);
            Row(index, range, result);
            return result;
        }

        /// <summary>
        /// Copies the given <see cref="Range"/> of a row into an <see cref="Vector"/>.
        /// </summary>
        /// <param name="index">The column to copy.</param>
        /// <param name="range">The <see cref="Range"/> of elements to copy.</param>
        /// <param name="result">The <see cref="Vector"/> to copy the row into.</param>
        /// <returns>A <see cref="Vector"/> containing the copied elements.</returns>
        /// <exception cref="ArgumentNullException">if the result matrix is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <b>result.Count != (range.End - range.Start) / range.Stride + 1</b>.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="range"/> specifies a range
        /// outside the row's range. </exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="index"/> is negative,
        /// or greater than or equal to the number of rows.</exception>    
        public override void Row(int index, Range range, Vector result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }
            if (index >= Rows
                || index < 0)
            {
                throw new ArgumentOutOfRangeException("index", Strings.OutOfRangeException);
            }
            if (range.Start >= Columns
                || range.End >= Columns)
            {
                throw new ArgumentOutOfRangeException("range", Strings.OutOfRangeException);
            }
            
            // Clear the result vector
            result.Clear();

            for (int i = range.Start, k = 0; i <= range.End; i += range.Stride, k++)
            {
                int itemIndex;
                if (!FindItem(index, i, out itemIndex))
                {
                    continue;
                }

                result[k] = mValues[itemIndex];
            }
        }

        /// <summary>
        /// Copies the elements of this matrix to the given matrix.
        /// </summary>
        /// <param name="target">The matrix to copy values into.</param>
        /// <exception cref="ArgumentNullException">if target is null.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if target is not the same size as this matrix.</exception>
        public override void CopyTo(Matrix target)
        {
            if (target == null)
            {
                throw new ArgumentNullException("target", Strings.NullParameterException);
            }

            if ((target.Rows != Rows) || (target.Columns != Columns))
            {
                throw new NotConformableException("target", Strings.ParameterNotConformable);
            }
            
            // Clear the target matrix
            target.Clear();

            SparseMatrix matrix = target as SparseMatrix;

            if (matrix != null)
            {
                // The matrix is a sparse matrix. Do a quick block copy of the
                // data.
                CopyToSparseMatrix(matrix);
            }
            else
            {
                // Do a normal copy
                // iterate over the rows
                for (int i = 0; i < mRowIndices.Length; i++)
                {
                    // Get the begin / end index for the current row
                    int begin;
                    int end;
                    FindRowStartAndEnd(i, out begin, out end);

                    // Get the values for the current row
                    if (begin == end)
                    {
                        // Begin and end are equal. There are no values in the row
                        // Move to the next row
                        continue;
                    }

                    for (int j = begin; j < end; j++)
                    {
                        // We could use the ValueAt method.
                        // target.ValueAt(i, mColumnIndices[j], mValues[j]);

                        target[i, mColumnIndices[j]] = mValues[j];
                    }
                }
            }
        }

        private void CopyToSparseMatrix(SparseMatrix matrix)
        {
            // copy the row indices
            Buffer.BlockCopy(mRowIndices, 0, matrix.mRowIndices, 0, mRowIndices.Length * Constants.SizeOfInt);

            // copy the colum indices. Only copy the actual data items. Note that the array
            // may be bigger than the number of items we copy.
            matrix.mColumnIndices = new int[mValueCount];
            Buffer.BlockCopy(mColumnIndices, 0, matrix.mColumnIndices, 0, mValueCount * Constants.SizeOfInt);

            // copy the data. Only copy the actual data items. Note that the array
            // may be bigger than the number of items we copy.
            matrix.mValueCount = mValueCount;
            matrix.mValues = new double[mValueCount];
            Buffer.BlockCopy(mValues, 0, matrix.mValues, 0, mValueCount * Constants.SizeOfDouble);
        }

        /// <summary>
        /// Copies the transpose of this matrix into the result matrix.
        /// </summary>
        /// <param name="result">the result of the transpose.</param>
        /// <exception cref="ArgumentNullException">if the result matrix is null.</exception> 
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the result matrix's dimensions are not this.Columns x this.Rows.</exception>
        public override void Transpose(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if ((result.Rows != Columns) || (result.Columns != Rows))
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            
            // Clear the result matrix
            result.Clear();

            // Do an 'inverse' CopyTo
            // iterate over the rows
            for (int i = 0; i < mRowIndices.Length; i++)
            {
                // Get the begin / end index for the current row
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);

                // Get the values for the current row
                if (begin == end)
                {
                    // Begin and end are equal. There are no values in the row
                    // Move to the next row
                    continue;
                }

                for (int j = begin; j < end; j++)
                {
                    // We could use the ValueAt method.
                    // target.ValueAt(i, mColumnIndices[j], mValues[j]);

                    result[mColumnIndices[j], i] = mValues[j];
                }
            }
        }

        /// <summary>
        /// Copies the values of the given <see cref="Vector"/> to the diagonal.
        /// </summary>
        /// <param name="source">The vector to copy the values from. The length of the vector should be
        /// Min(Rows, Columns).</param>
        /// <exception cref="ArgumentNullException">If <paramref name="source"/> is <b>null</b>.</exception>   
        /// <exception cref="NotConformableException">If the length of <paramref name="source"/> does not
        /// equal Min(Rows, Columns).</exception>
        /// <remarks>For non-square matrices, the elements of <paramref name="source"/> are copied to
        /// this[i,i].</remarks>       
        public override void SetDiagonal(Vector source)
        {
            if (source == null)
            {
                throw new ArgumentNullException("source", Strings.NullParameterException);
            }

            int min = System.Math.Min(Rows, Columns);
            if (source.Count != min)
            {
                throw new NotConformableException("source", Strings.ParameterNotConformable);
            }

            for (int i = 0; i < source.Count; i++)
            {
                int index;
                if (FindItem(i, i, out index))
                {
                    if (System.Math.Abs(source[i]) < Double.Epsilon)
                    {
                        DeleteItem(i, i);
                    }
                    else
                    {
                        mValues[index] = source[i];
                    }
                }
                else
                {
                    InsertValue(i, i, source[i]);
                }
            }
        }

        /// <summary>
        /// Copies the values of the given array to the diagonal.
        /// </summary>
        /// <param name="source">The array to copy the values from. The length of the array should be
        /// Min(Rows, Columns).</param>
        /// <exception cref="ArgumentNullException">if <paramref name="source"/> is <b>null</b>.</exception>        
        /// <exception cref="NotConformableException">If the length of <paramref name="source"/> does not
        /// equal Min(Rows, Columns).</exception>
        /// <remarks>For non-square matrices, the elements of <paramref name="source"/> are copied to
        /// this[i,i].</remarks>    
        public override void SetDiagonal(double[] source)
        {
            if (source == null)
            {
                throw new ArgumentNullException("source", Strings.NullParameterException);
            }

            int min = System.Math.Min(Rows, Columns);
            if (source.Length != min)
            {
                throw new NotConformableException("source", Strings.ParameterNotConformable);
            }

            for (int i = 0; i < source.Length; i++)
            {
                int index;
                if (FindItem(i, i, out index))
                {
                    if (System.Math.Abs(source[i]) < Double.Epsilon)
                    {
                        DeleteItem(i, i);
                    }
                    else
                    {
                        mValues[index] = source[i];
                    }
                }
                else
                {
                    InsertValue(i, i, source[i]);
                }
            }
        }

        /// <summary>
        /// Puts the the lower triangle of this matrix into the result matrix.
        /// </summary>
        /// <param name="result">where to store the the lower triangle.</param>
        /// <exception cref="ArgumentNullException">if the result matrix is null.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the result matrix's dimensions are not the same as this matrix.</exception>
        public override void LowerTriangle(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if ((result.Rows != Rows)
                || (result.Columns != Columns))
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            // Clear the result matrix
            result.Clear();
            
            for (int i = 0; i < mRowIndices.Length; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);

                for (int j = begin; j < end; j++)
                {
                    if (i >= mColumnIndices[j])
                    {
                        result[i, mColumnIndices[j]] = mValues[j];
                    }
                }
            }
        }

        /// <summary>
        /// Puts the the lower triangle of this matrix into the result matrix. The result matrix doesn't 
        /// contain the diagonal elements of this matrix.
        /// </summary>
        /// <param name="result">where to store the the lower triangle.</param>
        /// <exception cref="ArgumentNullException">if the result matrix is null.</exception>  
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the result matrix's dimensions are not the same as this matrix.</exception>
        public override void StrictlyLowerTriangle(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if ((result.Rows != Rows)
                || (result.Columns != Columns))
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            // Clear the result matrix
            result.Clear();
            
            for (int i = 0; i < mRowIndices.Length; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);

                for (int j = begin; j < end; j++)
                {
                    if (i > mColumnIndices[j])
                    {
                        result[i, mColumnIndices[j]] = mValues[j];
                    }
                }
            }
        }

        /// <summary>
        /// Puts the the upper triangle of this matrix into the result matrix.
        /// </summary>
        /// <param name="result">where to store the the upper triangle.</param>
        /// <exception cref="ArgumentNullException">if the result matrix is null.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the result matrix's dimensions are not the same as this matrix.</exception>
        public override void UpperTriangle(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if ((result.Rows != Rows)
                || (result.Columns != Columns))
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            // Clear the result matrix
            result.Clear();
            
            for (int i = 0; i < mRowIndices.Length; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);

                for (int j = begin; j < end; j++)
                {
                    if (i <= mColumnIndices[j])
                    {
                        result[i, mColumnIndices[j]] = mValues[j];
                    }
                }
            }
        }

        /// <summary>
        /// Puts the the upper triangle of this matrix into the result matrix. The result matrix doesn't 
        /// contain the diagonal elements of this matrix.
        /// </summary>
        /// <param name="result">where to store the the upper triangle.</param>
        /// <exception cref="ArgumentNullException">if the result matrix is null.</exception>          
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the result matrix's dimensions are not the same as this matrix.</exception>
        public override void StrictlyUpperTriangle(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if ((result.Rows != Rows)
                || (result.Columns != Columns))
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            
            // Clear the result matrix
            result.Clear();
            
            for (int i = 0; i < mRowIndices.Length; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);

                for (int j = begin; j < end; j++)
                {
                    if (i < mColumnIndices[j])
                    {
                        result[i, mColumnIndices[j]] = mValues[j];
                    }
                }
            }
        }

        /// <summary>
        /// Copies a sub-matrix of this matrix into a given matrix.
        /// </summary>
        /// <param name="rowRange">The <see cref="Range"/> of rows to copy into the given matrix.</param>
        /// <param name="columnRange">The <see cref="Range"/> of columns to copy into the given matrix.</param>
        /// <param name="result">The matrix to copy the sub-matrix into.</param>
        /// <exception cref="ArgumentNullException">if the result matrix is null.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the result matrix's dimensions are not the same as those defined by the <b>Ranges</b>.</exception>
        public override void SubMatrix(Range rowRange, Range columnRange, Matrix result)
        {
            if (rowRange.Start >= Rows
                || rowRange.End >= Rows)
            {
                throw new ArgumentOutOfRangeException("rowRange", Strings.OutOfRangeException);
            }

            if (columnRange.Start >= Columns
                || columnRange.End >= Columns)
            {
                throw new ArgumentOutOfRangeException("columnRange", Strings.OutOfRangeException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }
            
            // Clear the result matrix
            result.Clear();

            int m = (rowRange.End - rowRange.Start) / rowRange.Stride + 1;
            int n = (columnRange.End - columnRange.Start) / columnRange.Stride + 1;
            if ((result.Rows != m) || (result.Columns != n))
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            for (int i = rowRange.Start, row = 0; i <= rowRange.End; i += rowRange.Stride, row++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    // check if the column index is in the range
                    if ((((mColumnIndices[j] - columnRange.Start) % columnRange.Stride) == 0) &&
                        (mColumnIndices[j] <= columnRange.End))
                    {
                        int column = ((mColumnIndices[j] - columnRange.Start) / columnRange.Stride);
                        result[row, column] = mValues[j];
                    }
                }
            }
        }

        /// <summary>
        /// Copies the value of a <b>Matrix</b> into this <b>Matrix</b> stating at the given
        /// <paramref name="row"/> and <paramref name="column"/>.
        /// </summary>
        /// <param name="row">The row of this <b>Matrix</b> to start copying the values to.</param>
        /// <param name="column">The column of this <b>Matrix</b> to start copying the values to.</param>
        /// <param name="subMatrix">The <b>Matrix</b> to copy the values from.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="subMatrix"/> is null.</exception>
        /// <exception cref="ArgumentOutOfRangeException">
        /// <list type="bullet">
        /// <item><term>row</term><description>The given row is greater than or equal to the number of rows of this <b>Matrix</b>.</description></item>
        /// <item><term>column</term><description>The given column is greater than or equal to the number of columns of this <b>Matrix</b>.</description></item>
        /// <item><term>subMatrix</term><description>The sub-matrix does not fit into this <b>Matrix</b>.</description></item>
        /// </list></exception>
        public override void SetSubMatrix(int row, int column, Matrix subMatrix)
        {
            if (subMatrix == null)
            {
                throw new ArgumentNullException("subMatrix", Strings.NullParameterException);
            }

            if (row >= Rows)
            {
                throw new ArgumentOutOfRangeException("row", Strings.OutOfRangeException);
            }

            if (column >= Columns)
            {
                throw new ArgumentOutOfRangeException("column", Strings.OutOfRangeException);
            }

            if (row + subMatrix.Rows > Rows || column + subMatrix.Columns > Columns)
            {
                throw new ArgumentOutOfRangeException("subMatrix", Strings.OutOfRangeException);
            }

            // Want to do a rowwise insert but we can't. The main problem is that
            // the row begin & end may change if items are inserted or deleted.
            for (int i = 0; i < subMatrix.Rows; i++)
            {
                for (int j = 0; j < subMatrix.Columns; j++)
                {
                    int index;
                    if (FindItem(i + row, j + column, out index))
                    {
                        if (System.Math.Abs(subMatrix[i, j]) < Double.Epsilon)
                        {
                            DeleteItem(i + row, j + column);
                        }
                        else
                        {
                            mValues[index] = subMatrix[i, j];
                        }
                    }
                    else
                    {
                        InsertValue(i + row, j + column, subMatrix[i, j]);
                    }
                }
            }
        }

        /// <summary>
        /// Returns this matrix as a multidimensional array.
        /// </summary>
        /// <returns>a multidimensional containing the values of this matrix.</returns>        
        [SuppressMessage("Microsoft.Performance", "CA1814:PreferJaggedArraysOverMultidimensional", MessageId = "Return"), SuppressMessage("Microsoft.Performance", "CA1814:PreferJaggedArraysOverMultidimensional", MessageId = "Body")]
        public override double[,] ToArray()
        {
            double[,] result = new double[Rows, Columns];
            for (int rowIndex = 0; rowIndex < mRowIndices.Length; rowIndex++)
            {
                int columnStart;
                int columnEnd;
                FindRowStartAndEnd(rowIndex, out columnStart, out columnEnd);

                for (int columnIndex = columnStart; columnIndex < columnEnd; columnIndex++)
                {
                    result[rowIndex, mColumnIndices[columnIndex]] = mValues[columnIndex];
                }
            }
            return result;
        }

        /// <summary>
        /// Adds a scalar to each element in the matrix overwriting the values of this matrix.
        /// </summary>
        /// <param name="scalar">the scalar to add.</param>        
        public override void Add(double scalar)
        {
            // First check if a is equal to zero. If so
            // then we can just exit
            if (System.Math.Abs(scalar) < Double.Epsilon)
            {
                return;
            }

            for (int i = 0; i < Rows; i++)
            {
                for (int j = 0; j < Columns; j++)
                {
                    int index;
                    if (FindItem(i, j, out index))
                    {
                        // check if a is the negative of mValues[index].
                        // If so then we just delete the value at index
                        if (System.Math.Abs(mValues[index] + scalar) < Double.Epsilon)
                        {
                            DeleteItem(i, j);
                        }
                        else
                        {
                            mValues[index] += scalar;
                        }
                    }
                    else
                    {
                        InsertValue(i, j, scalar);
                    }
                }
            }
        }

        /// <summary>
        /// Adds a scalar to each element in the matrix and places the results into the result matrix.
        /// </summary>
        /// <param name="scalar">the scalar to add.</param>
        /// <param name="result">the results of the addition.</param>
        /// <exception cref="ArgumentNullException">if the result matrix is null.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the result matrix's dimensions are not the same as this matrix.</exception>
        public override void Add(double scalar, Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if ((result.Rows != Rows) || (result.Columns != Columns))
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            
            // Clear the result matrix
            result.Clear();

            for (int i = 0; i < Rows; i++)
            {
                for (int j = 0; j < Columns; j++)
                {
                    result[i, j] = this[i, j] + scalar;
                }
            }
        }

        /// <summary>
        /// Adds another matrix to this matrix overwriting the values of this matrix.
        /// </summary>
        /// <param name="other">the matrix to add to this matrix.</param>
        /// <exception cref="ArgumentNullException">if the other matrix is null.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the two matrices don't have the same dimensions.</exception>
        public override void Add(Matrix other)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other", Strings.NullParameterException);
            }

            if ((other.Rows != Rows) || (other.Columns != Columns))
            {
                throw new NotConformableException("other", Strings.ParameterNotConformable);
            }

            SparseMatrix matrix = other as SparseMatrix;
            if (matrix != null)
            {
                AddSparseMatrixToSelf(matrix);
            }
            else
            {
                for (int i = 0; i < other.Rows; i++)
                {
                    for (int j = 0; j < other.Columns; j++)
                    {
                        int index;
                        if (FindItem(i, j, out index))
                        {
                            // check if a is the negative of mValues[index].
                            // If so then we just delete the value at index
                            if (System.Math.Abs(mValues[index] + other[i, j]) < Double.Epsilon)
                            {
                                DeleteItem(i, j);
                            }
                            else
                            {
                                mValues[index] += other[i, j];
                            }
                        }
                        else
                        {
                            InsertValue(i, j, other[i, j]);
                        }
                    }
                }
            }
        }

        private void AddSparseMatrixToSelf(SparseMatrix matrix)
        {
            for (int i = 0; i < matrix.Rows; i++)
            {
                int begin;
                int end;
                matrix.FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    int index;
                    if (FindItem(i, matrix.mColumnIndices[j], out index))
                    {
                        // check if a is the negative of mValues[index].
                        // If so then we just delete the value at index
                        if (System.Math.Abs(mValues[index] + matrix.mValues[j]) < Double.Epsilon)
                        {
                            DeleteItem(i, matrix.mColumnIndices[j]);
                        }
                        else
                        {
                            mValues[index] += matrix.mValues[j];
                        }
                    }
                    else
                    {
                        InsertValue(i, matrix.mColumnIndices[j], matrix.mValues[j]);
                    }
                }
            }
        }

        /// <summary>
        /// Negates each element of this matrix.
        /// </summary>        
        public override void Negate()
        {
            for (int i = 0; i < mValueCount; i++)
            {
                mValues[i] = -mValues[i];
            }
        }

        /// <summary>
        /// Subtracts another matrix from this matrix overwriting the values of this matrix.
        /// </summary>
        /// <param name="other">the matrix to subtract.</param>
        /// <exception cref="ArgumentNullException">if the other matrix is null.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the two matrices don't have the same dimensions.</exception>
        public override void Subtract(Matrix other)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other", Strings.NullParameterException);
            }

            if ((other.Rows != Rows) || (other.Columns != Columns))
            {
                throw new NotConformableException("other", Strings.ParameterNotConformable);
            }

            SparseMatrix matrix = other as SparseMatrix;
            if (matrix != null)
            {
                SubstractSparseMatrixFromSelf(matrix);
            }
            else
            {
                for (int i = 0; i < other.Rows; i++)
                {
                    for (int j = 0; j < other.Columns; j++)
                    {
                        int index;
                        if (FindItem(i, j, out index))
                        {
                            // check if a is the positive of mValues[index].
                            // If so then we just delete the value at index
                            if (System.Math.Abs(mValues[index] - other[i, j]) < Double.Epsilon)
                            {
                                DeleteItem(i, j);
                            }
                            else
                            {
                                mValues[index] -= other[i, j];
                            }
                        }
                        else
                        {
                            InsertValue(i, j, -other[i, j]);
                        }
                    }
                }
            }
        }

        private void SubstractSparseMatrixFromSelf(SparseMatrix matrix)
        {
            for (int i = 0; i < matrix.Rows; i++)
            {
                int begin;
                int end;
                matrix.FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    int index;
                    if (FindItem(i, matrix.mColumnIndices[j], out index))
                    {
                        // check if a is the positive of mValues[index].
                        // If so then we just delete the value at index
                        if (System.Math.Abs(mValues[index] - matrix.mValues[j]) < Double.Epsilon)
                        {
                            DeleteItem(i, matrix.mColumnIndices[j]);
                        }
                        else
                        {
                            mValues[index] -= matrix.mValues[j];
                        }
                    }
                    else
                    {
                        InsertValue(i, matrix.mColumnIndices[j], -matrix.mValues[j]);
                    }
                }
            }
        }

        /// <summary>
        /// Multiplies each element of this matrix with a scalar overwriting the values of this matrix.
        /// </summary>
        /// <param name="scalar">the scalar to multiply with.</param>
        public override void Multiply(double scalar)
        {
            // if double = 0.0 then just reset the
            // value count to zero
            if (System.Math.Abs(scalar) < Double.Epsilon)
            {
                Clear();
                return;
            }

            for (int i = 0; i < mValueCount; i++)
            {
                mValues[i] *= scalar;
            }
        }

        /// <summary>
        /// Multiplies this matrix with another matrix and places the results into the result matrix.
        /// </summary>
        /// <param name="other">the matrix to multiply with.</param>
        /// <param name="result">the result of the multiplication.</param>
        /// <exception cref="ArgumentNullException">if the other matrix is null.</exception>
        /// <exception cref="ArgumentNullException">if the result matrix is null.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if this.Columns != other.Rows</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the result matrix's dimensions are not the this.Rows x other.Columns.</exception>
        public override void Multiply(Matrix other, Matrix result)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (Columns != other.Rows)
            {
                throw new NotConformableException("other", Strings.ParameterNotConformable);
            }

            if ((result.Rows != Rows) || (result.Columns != other.Columns))
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            
            // Clear the result matrix
            result.Clear();

            // The data is stored in row wise fashion.
            // So compute the values in row wise fashion
            DenseVector columnVector = (DenseVector)VectorBuilder.CreateVector(other.Rows, VectorType.Dense);

            for (int rowIndex = 0; rowIndex < Rows; rowIndex++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(rowIndex, out begin, out end);
                for (int columIndex = 0; columIndex < other.Columns; columIndex++)
                {
                    // clear the vector so that we start with an empty vector every time.
                    // Currently the SparseMatrix.Column doesn't clear the vector.
                    columnVector.Clear();
                    other.Column(columIndex, columnVector);

                    double sum = 0;
                    for (int k = begin; k < end; k++)
                    {
                        sum += mValues[k] * columnVector[mColumnIndices[k]];
                    }
                    result[rowIndex, columIndex] = sum;
                }
            }
        }

        /// <summary>
        /// Multiples this matrix with a vector and places the results into the result matrix.
        /// </summary>
        /// <param name="rightSide">the vector to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        /// <exception cref="ArgumentNullException">if rightSide is null.</exception>
        /// <exception cref="ArgumentNullException">if the result matrix is null.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the result vector size is not rightSide.Count.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the this.Columns != rightSide.Count.</exception>
        public override void Multiply(Vector rightSide, Vector result)
        {
            if (rightSide == null)
            {
                throw new ArgumentNullException("rightSide", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (rightSide.Count != Columns)
            {
                throw new NotConformableException("rightSide", Strings.ParameterNotConformable);
            }

            if (result.Count != Rows)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            
            // Clear the result vector
            result.Clear();

            // For each row in the matrix do:
            // 1) find out where the row starts and ends
            // 2) for each colum number find the appropriate value in a
            // 3) Multiply a and the value in the matrix
            for (int i = 0; i < Rows; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);

                result[i] = 0;
                for (int j = begin; j < end; j++)
                {
                    result[i] += mValues[j] * rightSide[mColumnIndices[j]];
                }
            }
        }

        /// <summary>
        /// Multiplies this matrix with a vector and places the results into the result matrix. Note the operation is vector * matrix.
        /// </summary>
        /// <param name="leftSide">The vector to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        /// <exception cref="ArgumentNullException">If leftSide is <b>null</b>.</exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <b>result.Count != this.Columns</b>.</exception>
        /// <exception cref="NotConformableException">If <b>this.Rows != leftSide.Count</b>.</exception>
        public override void LeftMultiply(Vector leftSide, Vector result)
        {
            if (leftSide == null)
            {
                throw new ArgumentNullException("rightSide", Strings.NullParameterException);
            }

            if (Rows != leftSide.Count)
            {
                throw new NotConformableException("leftSide", Strings.ParameterNotConformable);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (Columns != result.Count)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            
            // Clear the result vector
            result.Clear();

            foreach (KeyValuePair<int, Vector> column in ColumnEnumerator())
            {
                result[column.Key] = leftSide.DotProduct(column.Value);
            }
        }

        /// <summary>
        /// Stacks one this matrix on top of the given matrix, return in the combined matrix.
        /// </summary>
        /// <param name="lower">the matrix to stack this matrix upon.</param>
        /// <returns>the combined matrix.</returns>
        /// <exception cref="ArgumentNullException">if lower is null.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if upper.Columns != lower.Columns.</exception>
        public override Matrix Stack(Matrix lower)
        {
            if (lower == null)
            {
                throw new ArgumentNullException("lower", Strings.NullParameterException);
            }

            if (lower.Columns != Columns)
            {
                throw new NotConformableException("lower", Strings.ParameterNotConformable);
            }

            Matrix result = MatrixBuilder.CreateMatrix(Rows + lower.Rows, Columns, MatrixType.Sparse);
            Stack(lower, result);
            return result;
        }

        /// <summary>
        /// Stacks one this matrix on top of the given matrix and places the result into the result matrix.
        /// </summary>
        /// <param name="lower">the matrix to stack this matrix upon.</param>
        /// <param name="result">the combined matrix.</param>
        /// <exception cref="ArgumentNullException">if lower is null.</exception>
        /// <exception cref="ArgumentNullException">if the result matrix is null.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if upper.Columns != lower.Columns.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the result matrix's dimensions are not (this.Rows + lower.rows) x this.Columns.</exception>
        public override void Stack(Matrix lower, Matrix result)
        {
            if (lower == null)
            {
                throw new ArgumentNullException("lower", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (lower.Columns != Columns)
            {
                throw new NotConformableException("lower", Strings.ParameterNotConformable);
            }

            if ((result.Rows != (Rows + lower.Rows)) || (result.Columns != Columns))
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            
            // Clear the result matrix
            result.Clear();

            SparseMatrix lowerMatrix = lower as SparseMatrix;
            SparseMatrix resultMatrix = result as SparseMatrix;
            if ((lowerMatrix != null) && (resultMatrix != null))
            {
                StackSparseMatricesIntoNewSparseMatrix(lowerMatrix, resultMatrix);
            }
            else
            {
                if (lowerMatrix != null)
                {
                    StackSparseMatrices(lowerMatrix, result);
                }
                else
                {
                    for (int i = 0; i < Rows; i++)
                    {
                        int begin;
                        int end;
                        FindRowStartAndEnd(i, out begin, out end);
                        for (int j = begin; j < end; j++)
                        {
                            result[i, mColumnIndices[j]] = mValues[j];
                        }
                    }

                    for (int i = 0; i < lower.Rows; i++)
                    {
                        for (int j = 0; j < lower.Columns; j++)
                        {
                            result[i + Rows, j] = lower[i, j];
                        }
                    }
                }
            }
        }

        private void StackSparseMatricesIntoNewSparseMatrix(SparseMatrix lower, SparseMatrix result)
        {
            // Have all sparse matrices. Just concenate the values array's
            result.mValues = new double[mValueCount + lower.mValueCount];
            result.mValueCount = mValueCount + lower.mValueCount;
            result.mColumnIndices = new int[mValueCount + lower.mValueCount];


            // Copy the row data
            Buffer.BlockCopy(mRowIndices, 0, result.mRowIndices, 0, mRowIndices.Length * Constants.SizeOfInt);
            // Now add the row data from the lower matrix. We can't copy this because we have to change the
            // values.
            for (int i = 0; i < lower.mRowIndices.Length; i++)
            {
                result.mRowIndices[i + mRowIndices.Length] = mValueCount + lower.mRowIndices[i];
            }

            // Copy the column data
            Buffer.BlockCopy(mColumnIndices, 0, result.mColumnIndices, 0, mValueCount * Constants.SizeOfInt);
            Buffer.BlockCopy(lower.mColumnIndices, 0, result.mColumnIndices, mValueCount * Constants.SizeOfInt,
                             lower.mValueCount * Constants.SizeOfInt);

            // Copy the values
            Buffer.BlockCopy(mValues, 0, result.mValues, 0, mValueCount * Constants.SizeOfDouble);
            Buffer.BlockCopy(lower.mValues, 0, result.mValues, mValueCount * Constants.SizeOfDouble, lower.mValueCount * Constants.SizeOfDouble);
        }

        private void StackSparseMatrices(SparseMatrix lower, Matrix result)
        {
            for (int i = 0; i < Rows; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result[i, mColumnIndices[j]] = mValues[j];
                }
            }

            for (int i = 0; i < lower.Rows; i++)
            {
                int begin;
                int end;
                lower.FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result[i + Rows, lower.mColumnIndices[j]] = lower.mValues[j];
                }
            }
        }

        /// <summary>
        /// Concatenates this matrix with the given matrix and places the result into the result matrix.
        /// </summary>
        /// <param name="right">The matrix to concatenate.</param>
        /// <param name="result">The combined matrix.</param>
        /// <exception cref="ArgumentNullException">If right is <b>null</b>.</exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <b>this.Rows != right.Rows</b>.</exception>
        /// <exception cref="NotConformableException">If the result matrix's dimensions are not this.Rows x (this.Columns + right.Columns).</exception>
        public override void Append(Matrix right, Matrix result)
        {
            if (right == null)
            {
                throw new ArgumentNullException("right", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (right.Rows != Rows)
            {
                throw new NotConformableException("right", Strings.ParameterNotConformable);
            }

            if ((result.Rows != Rows) || (result.Columns != (Columns + right.Columns)))
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            
            // Clear the result matrix
            result.Clear();

            SparseMatrix lowerMatrix = right as SparseMatrix;
            SparseMatrix resultMatrix = result as SparseMatrix;
            if ((lowerMatrix != null) && (resultMatrix != null))
            {
                AppendSparseMatricesIntoNewSparseMatrix(lowerMatrix, resultMatrix);
            }
            else
            {
                if (lowerMatrix != null)
                {
                    AppendSparseMatrices(lowerMatrix, result);
                }
                else
                {
                    for (int i = 0; i < Rows; i++)
                    {
                        int begin;
                        int end;
                        FindRowStartAndEnd(i, out begin, out end);
                        for (int j = begin; j < end; j++)
                        {
                            result[i, mColumnIndices[j]] = mValues[j];
                        }
                    }

                    for (int i = 0; i < right.Rows; i++)
                    {
                        for (int j = 0; j < right.Columns; j++)
                        {
                            result[i, j + Columns] = right[i, j];
                        }
                    }
                }
            }
        }

        private void AppendSparseMatricesIntoNewSparseMatrix(SparseMatrix right, SparseMatrix result)
        {
            // Have all sparse matrices. Just concenate the values array's
            result.mValues = new double[mValueCount + right.mValueCount];
            result.mValueCount = mValueCount + right.mValueCount;
            result.mColumnIndices = new int[mValueCount + right.mValueCount];

            int rowStart = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                int begin;
                int end;

                FindRowStartAndEnd(i, out begin, out end);
                int length = 0;
                for (int j = begin; j < end; j++)
                {
                    result.mValues[rowStart + length] = mValues[j];
                    result.mColumnIndices[rowStart + length] = mColumnIndices[j];
                    length++;
                }

                right.FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result.mValues[rowStart + length] = right.mValues[j];
                    result.mColumnIndices[rowStart + length] = right.mColumnIndices[j] + Columns;
                    length++;
                }

                result.mRowIndices[i] = rowStart;
                rowStart += length;
            }
        }

        private void AppendSparseMatrices(SparseMatrix right, Matrix result)
        {
            for (int i = 0; i < Rows; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result[i, mColumnIndices[j]] = mValues[j];
                }
            }

            for (int i = 0; i < right.Rows; i++)
            {
                int begin;
                int end;
                right.FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result[i, Columns + right.mColumnIndices[j]] = right.mValues[j];
                }
            }
        }

        /// <summary>
        /// Diagonally stacks his matrix on top of the given matrix. The new matrix is a MxN matrix, 
        /// where M = this.Rows + lower.Rows and N = this.Columns + lower.Columns.
        /// The values of off the off diagonal matrices/blocks are set to zero.
        /// </summary>
        /// <param name="lower">the lower, right matrix.</param>
        /// <exception cref="ArgumentNullException">if lower is null.</exception>
        /// <returns>the combined matrix</returns>
        public override Matrix DiagonalStack(Matrix lower)
        {
            if (lower == null)
            {
                throw new ArgumentNullException("lower", Strings.NullParameterException);
            }

            Matrix result = MatrixBuilder.CreateMatrix(Rows + lower.Rows, Columns + lower.Columns, MatrixType.Sparse);
            DiagonalStack(lower, result);
            return result;
        }

        /// <summary>
        /// Diagonally stacks his matrix on top of the given matrix and places the combined matrix into the result matrix.
        /// </summary>
        /// <param name="lower">the lower, right matrix.</param>
        /// <param name="result">the combined matrix</param>
        /// <exception cref="ArgumentNullException">if lower is null.</exception>
        /// <exception cref="ArgumentNullException">if the result matrix is null.</exception>
        /// <exception cref="dnAnalytics.Exceptions.NotConformableException">if the result matrix's dimensions are not (this.Rows + lower.rows) x (this.Columns + lower.Columns).</exception>
        public override void DiagonalStack(Matrix lower, Matrix result)
        {
            if (lower == null)
            {
                throw new ArgumentNullException("lower", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if ((result.Rows != (Rows + lower.Rows)) || (result.Columns != (Columns + lower.Columns)))
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            
            // Clear the result matrix
            result.Clear();

            SparseMatrix lowerMatrix = lower as SparseMatrix;
            SparseMatrix resultMatrix = result as SparseMatrix;
            if ((lowerMatrix != null) && (resultMatrix != null))
            {
                DiagonalStackSparseMatricesIntoNewSparseMatrix(lowerMatrix, resultMatrix);
            }
            else
            {
                if (lowerMatrix != null)
                {
                    DiagonalStackSparseMatrices(lowerMatrix, result);
                }
                else
                {
                    for (int i = 0; i < Rows; i++)
                    {
                        int begin;
                        int end;
                        FindRowStartAndEnd(i, out begin, out end);
                        for (int j = begin; j < end; j++)
                        {
                            result[i, mColumnIndices[j]] = mValues[j];
                        }
                    }

                    for (int i = 0; i < lower.Rows; i++)
                    {
                        for (int j = 0; j < lower.Columns; j++)
                        {
                            result[i + Rows, j + Columns] = lower[i, j];
                        }
                    }
                }
            }
        }

        private void DiagonalStackSparseMatricesIntoNewSparseMatrix(SparseMatrix lower, SparseMatrix result)
        {
            // Have all sparse matrices. Just concenate the values array's
            result.mValues = new double[mValueCount + lower.mValueCount];
            result.mValueCount = mValueCount + lower.mValueCount;
            result.mColumnIndices = new int[mValueCount + lower.mValueCount];


            // Copy the row data
            Buffer.BlockCopy(mRowIndices, 0, result.mRowIndices, 0, mRowIndices.Length * Constants.SizeOfInt);
            // Now add the row data from the lower matrix. We can't copy this because we have to change the
            // values.
            for (int i = 0; i < lower.mRowIndices.Length; i++)
            {
                result.mRowIndices[i + mRowIndices.Length] = mValueCount + lower.mRowIndices[i];
            }

            // Copy the column data
            Buffer.BlockCopy(mColumnIndices, 0, result.mColumnIndices, 0, mValueCount * Constants.SizeOfInt);
            // Now add the row data from the lower matrix. We can't copy this because we have to change the
            // values.
            for (int i = 0; i < lower.mValueCount; i++)
            {
                result.mColumnIndices[i + mValueCount] = Columns + lower.mColumnIndices[i];
            }

            // Copy the values
            Buffer.BlockCopy(mValues, 0, result.mValues, 0, mValueCount * Constants.SizeOfDouble);
            Buffer.BlockCopy(lower.mValues, 0, result.mValues, mValueCount * Constants.SizeOfDouble, lower.mValueCount * Constants.SizeOfDouble);
        }

        private void DiagonalStackSparseMatrices(SparseMatrix lower, Matrix result)
        {
            for (int i = 0; i < Rows; i++)
            {
                int begin;
                int end;
                FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result[i, mColumnIndices[j]] = mValues[j];
                }
            }

            for (int i = 0; i < lower.Rows; i++)
            {
                int begin;
                int end;
                lower.FindRowStartAndEnd(i, out begin, out end);
                for (int j = begin; j < end; j++)
                {
                    result[i + Rows, Columns + lower.mColumnIndices[j]] = lower.mValues[j];
                }
            }
        }

        /// <summary>
        /// Returns a deep-copy clone of the <b>Matrix</b>.
        /// </summary>
        /// <returns>A deep-copy clone of the <b>Matrix</b>.</returns>
        public override Matrix Clone()
        {
            return new SparseMatrix(this);
        }
    }
}
